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Due to dwindling fossil fuel reserves, there is a global demand to increase 
the level of low-carbon based renewable energy resources (RERs) for 
electric power generation. Coupled with concerns that emissions from fossil 
fuels is leading to climate change with possible disastrous consequences, 
effort is seriously under way to discover the probable usage of RERs on 
large-scale without being integrated into an existing power grid. It is 
envisaged that such a large-scale RER power grid will operate solely on its 
own and expected to be stable and reliable comparable to a conventional 
power grid. The impact of wind power generation as part of the electric 
power grid is no longer negligible. Wind energy generation is one of the 
most established renewable energy resources to help ensure low carbon- 
based renewable energy (RE) self-sufficiency, yet it is also one of the most 
volatile RERs. Despite its disadvantages, wind power generation is expected 
to continue its strong growth in the coming years as result of high interest in 
clean energy to curb the global warming. Various studies are looking at the 
prospects for solely RER power grids for usage on large-scale. However, the 
issue has been the stability and the reliability of such power grids. 
Variability of power output, intermittency and load mismatch are wind 
farms’ unique characteristics potentially harmful to grid voltage stability. In 
response to this problem, A large-scale wind power system was modelled 
using stochastic approach and the results analyzed using Lyapunov method, 
matrix exponential and Hurwitz criterion to ascertain the stability behavior 
of an entirely 100% RE grid being envisaged for the near future. 
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1. INTRODUCTION 


To achieve the year 2030 agenda for sustainable development (SD) as well as the Paris agreement 
on climate change, energy availability and affordability is key. Effort is being made by several countries 
around the world towards the search for sustainable and renewable energy (RE) alternatives to supplement 
their energy requirement due to factors such as the increasing demand for energy, the decline in fossil fuel 
reserves, CO2 reduction and global climate change [1]. The effects of the rise in global temperature is of great 
concern. Greenhouse gas (GHG) emissions which is as result of human activities is seriously affecting the 
Earth’s climate. It has been reported that the Earth’s climate has warmed up between the range of 0.8 °C to 
1.2 °C since 1882 and may reach 1.5 °C between the year 2030 and the year 2052 unless drastic measures are 
put in place. One of the most visible effects of global warming has been the rapid decrease in glaciation 
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experienced over the years. It is being speculated in certain quarters that the rising global temperature being 
experienced in modern times is due, in part, to atmospheric pollution arising from human activities, which is 
fuelling the Earth’s GHG effect. Although the GHG effect is essential to the well-being of mankind, 
excessive GHG levels above their natural norm, the corresponding additional warming could threaten the 
sustainability of the planet as a whole. The prediction is that is as global warming progresses, sea levels will 
rise by over 400 mm by the year 2080 due to the combined effects of thermal expansion of the oceans and 
melting of the polar ice. This will negatively impact livelihood of mankind, with an over 80 million people 
threatened with flooding in the low-lying territories [2], [3]. 

Electricity generation from RERs was said to have increased by over 8% in the year 2021, 
accounting for more than half of the increase in overall electricity production worldwide. RERs were 
expected to provide 30% of electricity generation worldwide by the end of the year 2021, their biggest share 
in the power generation mixes since the beginning of the Industrial Revolution. The largest contribution to 
that growth was expected to come from solar photovoltaic (PV) and wind. Electricity generation from wind is 
expected to grow by 275 terawatt-hours, or around 17%, as compared to that of the year 2020. Electricity 
generation from solar PV is expected to increase by 145 terawatt-hours, an increment of about 18% as 
compared to that of the year 2020. Their combined output was expected to have hit about 2800 terawatt- 
hours at the end of the year 2021 [3], [4]. A number of studies are looking at the prospects for solely 
renewable energy (RE) power grids, on large-scale. RERs are likely to account for about 95% of the net 
increase in global electric power supply by the year 2025 [5], [6]. The simplest way to meet these ambitions 
is to reduce the emissions from fossil fuels by deploying large-scale RERs [6]. 

The literature provided interesting expose regarding stochastic modelling of large-scale RERs 
integrated electric power grids. However, stochastic modelling of entirely RER electric power grid on large- 
scale was not found. The recent ones on RERs integrated large-scale electric power grids include: [7] 
established stochastic transient stability analysis method by introducing stochastic differential equations into 
the traditional transient stability analysis method. Utilizing the theory of probability and statistics, the Monte 
Carlo simulation was utilized to improve the numerical simulation method for transient stability evaluation. 
Using formulated impedance network model of large-scale renewable energy bases with frequency coupling, 
the s-domain nodal admittance matrix (SNAM) was formulated [8]. The stability analysis of a 10-node 
system containing 4 PMSG-based wind farms was carried out using the SNAM method. Kessywa et al. [9] 
conducted dynamic modelling and control study for the assessment of large-scale wind and solar photovoltaic 
integration into existing power system. The model was implemented in MATLAB-based transient stability 
analysis toolbox in order to analyses the dynamic response of the renewable energy sources. The developed 
model assumed that the converter is the only means through which the renewable energy resources interfaced 
the existing the power system. A novel method was proposed for a wind power output scene simulation [10]. 
A genetic algorithm (GA) K-means was used to divide the wind farm into clusters where K is the cluster 
number. The wind power output of each cluster was computed using the wind turbine model. Power output 
scenes were simulated based on Markov chain Monte Carlo (MCMC) method. A probabilistic analysis 
approach was used to investigate the impact of stochastic uncertainty of grid-connected wind generation on 
power system small-signal stability. The approach can compute the probabilistic density function (PDF) of 
critical eigenvalues of a large-scale multiple source wind integrated power system [11]. 

Model-based and measurement-based analyses are the two main approaches for determining an 
electric power system’s behavior under operation. Measurement-based analysis is a real-time monitoring 
approach for electric power systems. It differs from model-based analysis. It employs phasor measurement 
units (PMUs) which are high-speed sensors that measure voltage and current synchro-phasors of the electric 
power system with high accuracy to obtain the necessary information from the electric power system. 
Measurement-based analysis provides accurate data and the prediction of the stability behavior of an electric 
power system based on synchro-phasor technology. When model-based analysis approach is used to 
determine the stability of an electric power system, the electric power system is described using mathematical 
models. In situations where the effect of probabilistic uncertainty is considered, stochastic dynamical models 
play an important role in such analyses. Among the analysis approaches that can be applied to the 
mathematical model include Lyapunov method, eigenvalue analysis, and matrix exponential method and 
Routh-Hurwitz stability criterion [12], [13]. For this research, matrix Lyapunov function, matrix exponential 
and Hurwitz criterion were used to examine the stability behavior of a large-scale wind power system 
modelled using stochastic approach 

Wind energy is expanding steadily and significantly over the world as a result of the demand for 
more low-carbon RERs in the generation of electric power. One of the important technologies that has been 
around for a while, can be placed in utility-scale facilities with large installed capacity, and can compete 
competitively with traditional generation sources is wind power [14]. A by-product of solar energy is wind. 
Air moves as a result of atmospheric pressure gradients, which cause wind. From high-pressure areas to low- 
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pressure areas, wind blows. The higher the wind speed and the greater the amount of wind power that can be 
harvested from the wind by a wind turbine (WT), the greater the air pressure gradient. It is on record that the 
total solar power received by the earth is around 1.8x10!! MW. Only 2% (3.6x10° MW) of the solar energy is 
converted into wind energy, and approximately 35% of wind energy is dissipated within 1000 m of the 
earth’s surface. The available wind power that can be converted into other forms of energy is estimated to be 
around 1.26x10° MW which is about 20 times the rate of the present global energy consumption. Wind 
energy in principle could meet the entire energy needs of the world if harnessed effectively and efficiently 
[15]. The wind speed that is available has a significant impact on the power output of a WT. The height of 
the terrain, time of year, season, and location all affect the wind speed. However, it is a resource that exhibits 
significant variability at almost all-time scales, from the length of turbulent gusts to daily and seasonal cycles 
to long-term changes brought on by climatic conditions. At these time scales, various challenges, such as 
issues with stability, power quality, and reliability, arise [14]. The wind speed affects the WT's efficiency and 
safety. Three sets of wind speeds have been established by the International Electrotechnical Commission 
(IEC) for use in assessing the effectiveness and dependability of a WT. As indicated in Figure 1, the WT 
manufacturers typically offer these three speeds for any WT they produce [15]-[18]. A WT's performance is 
indicated by its power curve. The power output is zero in the initial zone where the wind speed is below a 
cut-in speed, or threshold minimum. There is a quick increase in power output in the second region, which is 
located between the cut-in and the rated speed. Until the cut-off speed is reached, the third zone produces a 
constant output (rated). Beyond this speed (region 4), the turbine is shut off to safeguard its internal parts from 
strong winds; as a result, no electricity is produced in this area. A power system comprising only RERs only 
should be possible to operate in the near future on large-scale due to the dwindling fossil fuel resources. 
However, the challenge to overcome in this regard is the stability and reliability of such future large-scale power 
systems. This research aims at developing a model to assess the stability of such future RERs power systems. 


Power 
output 
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0 Veutin Vrated Veut off 


— 
Wind speed 


Figure 1. Typical wind power curve of a wind turbine 


2. MODEL FORMULATION 
This section describes how the model was formulated. 


2.1. Wind power generation system modelling based on chapman-kolmogorov (C-K) equation 

In stochastic model, the knowledge of the state variables is subject to uncertainties due to some 
inherent randomness. This type of model is suitable for applications where the state information cannot be 
perfectly predictable. Markov processes are often used to model randomness and have no memory. The 
differential form of the Chapman-Kolmogorov equation is the Chemical Master Equation (CME) for jump 
processes or diffusive processes. The CME describes the evolution of the probability distribution of states of 
a continuous-time Markov process and is commonly used to describe stochastic physical, chemical, or 
biological systems. Jump processes are characterised by discontinuous motion; thus, there is a bounded and 
non-vanishing transition probability per unit time. The ME is a stochastic model and is derived from the 
fundamental assumption on the dynamic properties of the underlying stochastic process referred to as the 
Markov property [19]. As a result of wind speed intermittency or variations, the power generation into the 
power grid is stochastic or random in character impacting the stability and reliability of the power grid. The 
fluctuation of the renewable power generation can be regarded as random variables which are real-valued 
functions defined on a sample or state space with the assignment of possible probabilities to the possible 
values of the random variables [20]. Due to the large-scale nature of the renewable power grid being 
considered, continuous random variable is used which can be easily handled analytically as compared to 
discrete random variables. Again, for large-scale modelling, continuous random variables are preferred [21]. 
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The model formulation is based on the modified wind power curve given in Figure 2. Mid speed was 
introduced between the cut-in speed and the rated output speed for the purposes of the modelling for the 
realization of four transition states. For this study, the intermittency nature and the changing pattern of the 
wind speed is transformed into a transition diagram as given in Figure 3 having four states. Transitions 
between the different states occur from time to time. One can represent the transitions among the various 
states as random events that do occur at some rates. Each transition state is assigned a probability, p and in 
between the states is the transition rates represented by œ and ô. Based on Figure 3, the 12 possible transition 
states are as given by Table 1. A stochastic process is a process that evolves probabilistically through various 
states attained at various times from a well-defined initial state xo at to. The probability density function (pdf) 
of the process depends on their respective states and times at particular moment. 
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Figure 2. Wind power curve of a wind turbine with mid- Figure 3. State transition diagram for wind 


speed introduced turbine systems 


Table 1. Wind speed variation transition states 
S/No. Current state Next state Transition rate 


1 0 1 61 
2 0 2 Ws 
3 0 3 Wa 
4 1 0 Wi 
5 1 2 62 
6 1 3 56 
7 2 0 ôs 
8 2 1 w2 
9 2 2 63 
10 3 0 64 
11 3 1 We 
12 3 2 w3 


The probability that the process will be at state x, at time, t,; then progress to state x2, at time, tz 
then goes through all other subsequent configurations or states to be at state x,, at time tn, given that it started 
from state Xo at time tọ is given by (1): 


Pala = (Gn tn), n-1 tn-1) Xn-2 tn-2). e. (X1,t1)|(%o, to)) (1) 


The probability of the system being at state, x, at time t,, multiplied by the probability of it being at state, x2 
at time, t, given that it has been at state, x, at t,, and at state xg, at time, tọ respectively can be expressed 
as [22], [23]: 
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Pz (x2, tz), (%1,t1)|(%o, to)) = P (Œ ty) Io, to) )Prj2(%2, t2)|(%1, t1), (Xo to)) (2) 


If all such probabilities are known, the stochastic process is fully specified [22], [23]. The ME for 
jump processes is given by (3). 


PED = § wale oe’ t) — w(e'lx)p (x, t)) dx’ (3) 


Where: p(x, t) = probability of being in state x at time, t; w(x|x’) = transition probability of moving from 
state x’ to state x per unit time. 

One-step stochastic process is continuous-time Markov process whose range consists of the integer 
k, and whose transition probability per unit time, w,,,, permits only jumps between adjacent sites [24]. Thus, 


= mtk 
Wkm = TmOkm-1 + Onbignas ) 


Wem = 1- (Tk + Gx) (4) 


where, 7, = probability per unit time that, being in state k, a jump occurs to site k — 1, g,= probability per 
unit time that, being in state k, a jump occurs to state k + 1. 
One-step transition probabilities are also given by conditional probabilities [25]: 


Qun) = Q(Xna1 = J /Xn = Din = 0,1,2..... (5) 
One-step transition probabilities combined in matrix of one-step transition probabilities, Q is given by: 


foo qoi Foz *** 
dio 911 12 | 
P= : : Eo cigs (6) 


fio qi qiz” 


where, Q;j= probability of a transition from state i to state j in one step or one time unit or in one jump. 

A wind turbine (WT) transitioning from one state to another can be regarded as a jump process or 
one-step transition. For this modelling, every transition state is taken to be positive since the WT is not 
supposed to be stationary hence, a current state depends on all other previous states which can be treated as a 
union of disjoint intervals based on the theory and concept of continuous random variable. For a continuous 
random variable X with probability density function (PDF) fx(x), the following holds [20]: 


fx(x) = 0 forall x € R T) 
Sof du = 1 (8) 
P(a < X < b) = Fy(b) — Fy(a) = J; fx (u)du (9) 
More generally, for a set A, P(X € A) = J, fy (w)du (10) 


To determine the probability of a continuous random variable X, for instance, P(X € [1, 2] U [5, 6]; 
one can write [20]: 


P(X € [1,2] U [5,6] = ff fr(wdut fo fę(u)du (11) 


Based in (3), the ME for jump processes is modified and applied to wind turbine (WT) power generation 
system which is given by (12): 


PED = Swale oe’ t) +wp, t)) dx’ (12) 


Assuming the system is at state x (state 1 or initial state) at time, t moving to state x — 1 (state 2 or 
next state) at time, t + At. The probability of moving from state x (state | or initial state) to state x — 1 (state 
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2 or next state) is given by QAt. The transition must occur during the time, At. For the transition to occur: 
i) The WT power generation system must be at state x (state 1 or initial state) with the probability, p,(t); and 
ii) Transition must occur from state x (state 1 or initial state) to state x — 1 (state 2 or next state) with the 
probability, QAt, where 2 = transition rate. 

The probability of transitioning from state x (state 1 or initial state) to state x — 1 (state 2 or next 
state) between the times ¢, t + At is the product of the two probabilities given by px(t)QAt. During that 
small time interval, At, only one change in state must occur, either x — 1 (decrement in wind speed) or x + 1 
(increment in wind speed). Transitioning to the next state depends on the previous state and the transition 
rate. Based on the above, dynamical equations can be developed for state x (state 1 or initial state), state x- 1 
(state 2 or next state) or x + 1 state (another possible state 2 or next state) as follows based on the rules of 
conditional probability [20], [22]. This can be translated into probability terms as follows: 


Poy (t + At) = Poy Ct) + P E)NA (13) 
That of state (x — 1) transition is given by: 

Poe-1y(t + At) = Poy Ct) + Pas) (()QAt (14) 
And that of state (x + 1) is also be given by: 

P+) t + At) = P+) (t) + P+ E)NA (15) 


Dividing through in (13), (14) and (15) respectively by At and taking the limit as At — 0, in (13), (14) and 
(15) becomes: 


we 2D, (16) 
dPp(x- 
ER fip) (17) 
aD (x 

amn = 2D (x41) (18) 


Generalizing and applying the principle to the entire four states with twelve transition diagrams of 
Figure 3 (with the transitions given in Table 1) which represents the movement pattern of WTs due to 
fluctuations in wind speed with 1, 2, 3, ............ „n number of wind farms comprising a number of WTs; the 
dynamical ME in stochastic transitional probability matrix-vector form: 


Py -( +O, +0) Oi 55) bu Pi 
d Pa = On -(@, +6, +651) Oy Ogi Pa 
dt Ps Os, ô» — (55, tO +ô) O; Ps, 
Pa Oy 551 Oy, (G4, + Og + O) Pa 
Pn -6y, +O +0) Oz ô; by Py 
d Py | Ôn - (Op +5, +6) On Op Py 
dt Px i Os dy — (55) +0 +559) Oy : 
Px Oy ôg Ox - (ô, Ẹ Og + O) Py 
Pin —(6,, +5, +,,) On 6, Oan Pin 
d| Panl On —(@,,+6,, +55,) Os, On Poy (19) 
dt Dn z Os, >, —(6;, +), +6;,) O, | Don 
Pan Osn Osn os, —(6,, + Ogn + On) Pan 
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In (19) is modified accordingly based in (12): 


Pu (6,,+@.,+@,,) Oi Ô Oy Pu 
d| Py |_ Ô (@ +ó, +8) Oy, Og, Pa 
dt| P Os, >, (6;,+@,,+03,) Os, Di, 

Dy Ox, Os 3, (Ô + Og + 0) Pa 

Po (a tOn + Oy) O Os) Oy Po 
d Pr» |_ Ov (O +04, +669) Oy, Op Po 
dt| Py Ory Oy (05) + @yy +z) Oxy Ps 

Pw Op Os by) (p + Oy + Dy) Px 

Pin (6, +s, + Osn) On 55 Sn Pin 
d Pon _ Ôn (,, +n + 56n) Oy, On Pon (20) 
dt| D;, Os, On, (6,,+@,, +03,) O, Ds, 

Pan Oin Ögn Ò, (6,, + On + Qs, ) Pan 


Where the first subscript represents the transition rate number, and the second subscript also represents the 
wind farm number. The state stochastic transition probability matrix for the first wind farm is therefore given 
by: 


(Ô +@;,+@,) Or, Os On 
a On (@,,+05,+6¢;) Oy, Oy (21) 
Os, >, (ô; +@,, +) O;, 
Or Ov Ô; (4, + Og + 0) 


The transition rate for the wind farms, which varies between 0.1 to 1.0, the least being 0.1 and the maximum 
being 1.0. The maximum being 1.0 is assumed for the wind farm model. In (21) results in stochastic 
transition probability matrix, U: 


(22) 


= = =. WwW 


1 
3 
1 
1 


SB uvu =e =. 
wv ee Re 


Even though large-scale WT power generation system with its associated size of state matrix, [U] should be 
in order of thousand or more, depending on the number of wind farms being considered that goes with its 
many eigenvalues, it is not necessary to compute all these eigenvalues [26]. 


2.2. Model assumptions 
The following assumptions were made during the development of the model: 
— The wind speed variation is a function of the weather conditions only; 
— A transition to the next state may occur at any point in time; 
— Only one transition can occur at a time (1-step transition); 
— The transition to another state is gradual over time (not immediate); 
— A transition from one state to the next state for all WTs constituting a wind farm occurs at the same time 
— The transition rates vary between 0.1 — 1.0; 
— The wind power plants are spread geographically to deal with the issue of short-term variability; 
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— The WTs operate in variable-speed mode, and the control system regulates the rotor speed to obtain peak 
efficiency in fluctuating winds; 

— The wind farms are sited in areas with high average wind speeds; and 

— The electrical load model is based on Markov property (current state does not depend on previous state). 


2.3. Determination of characteristic equation 
The stochastic transitional probability matrix, U given by (22): 


3 1 1 1 
1 3 1 1 22 
on (22) 
1 1 3 1 
1 1 1 3 
The matrix, U has non-trivial solution given by (23) which says that the values of à which satisfies: 
det|U — AI = 0| (23) 


Where: A = unknown eigenvalues of the matrix, U; I =n x n identity matrix. In general, for any n X n matrix U, 


my-A u2 Win 
eee. ae un (24) 
un] U yd) eeen ünn a 


Det U-Al) = (u m7” det}... a wae |E 


Based in (24) and (25), (22) is solved as follows for the eigenvalues or characteristic roots of the wind farm 1 
(WF-1) made up of a number of wind turbines as follows: 


ao 00 3 1 1 1 
0240 0 1 3 1 1 6 (26) 
0 020 1 1 3 1 
00 04 1 1 1 3 
3-A 1 1 1 
3-A 1 1 (27) 
= 0 
1 1 3-A 1 
1 1 1 3-A 
The determinant in (27) is given by (28): 
1.02°-18.045-128.044—46413 +912.0A7-928.01+384.0 (28) 
1.0A?—6.0A+8.0 
Therefore, the characteristic equations of wind farms 1, 2.......... n is given by: 
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Aa —124? +4847 —804,+48 =0 
A,* -124,3 +4847 —804, +48 = 0 


4,4 -124,3 +484,7 —804, +48 =0 (29) 


2.4. Lyapunov stability analysis approach 

To assess the stability of a particular system using energy function approach, the Lyapunov function 
and its derivative with respect to time is determined. Lyapunov function may include physical variables or 
the whole state variables of the system under consideration. The quadratic Lyapunov function and its 
derivative with respect to time are defined as follows for any linear system [28], [29]: 


x=xA (30) 
Its corresponding Lyapunov function is defined in matrix form as: 
V(x) = x? Px (31) 


Where, P is a real symmetric matrix. Its derivative is given by: 


— = = x (AP + PA 
de Goes (32) 
= x"(A’P + PA)x = x" Qx 


Where, x and x’ are the state variables’ vector and its transpose vector respectively. For any positive definite 
Hermitian or symmetric matrix Q, a positive definite Hermitian or symmetric matrix P, satisfying the 
following Lyapunov matrix equation given by (33) based on (32) [28]-[32]. 


A™P+PA= -Q (33) 


Where, A, P, Q = Q" ER; holds for some positive definite Q = Q7 > 0 and P = P7 > 0, then matrix A is stable. 
Where, Q” and P” are the transpose of the matrices Q and P respectively. 

To evaluate the stability of any dynamical system using the Lyapunov-based approach, the first step 
is to choose an appropriate Q matrix, which should be positive definite. If both P and Q are positive definite, 
then the system works in the stable mode. By defining the P, Q will be determined based on the system state 
matrix. Alternatively, by defining Q, P can be determined. For this particular study, Q is defined in a specific 
positive definite manner, then system stability is determined based on the analyzing the matrix, P. The 
quadratic function, V(x) based on matrix given by (22) is defined as [27]—[29], [33]: 


Pa Py Po Pua ||% 
Pn Pn Po Pu ||% 
P13 Po P3 Po ||% 
Pig Pa P34 Pu J %4 


V(x) =x Px = [x X, x; x4] 


(34) 


2 2 
V(x) = PX +2 DyX%Xy +2 pX% +2 DyXXy + Py Xp +2 PX% +2 Pa XX4 


2 2 
+ P33X3 +2 P34X3X4 + Puy X4 


Where, x is real vector and P is real symmetric matrix. Based in (33), Q is chosen as positive-definite unit 
matrix given by: 
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1 0 0 0 
01 0 0 (33) 
Q= 
0 0 1 0 
0 0 0 1 
The real symmetric n x n matrix P, is also given by: 
Py Pu Pi Pu 
Pi2 Px P23 Pa (36) 


Pis P23 P33 P3 
Pi4 Pa P34 Pag 


Where the matrix P’s elements are P11, Piz. . .P44. Substituting in (35) and (36) into (33) to obtain the matrix, 
P expressed by (36): 


—0.2083 0.0417 0.0417 0.0417 

0.0417 —0.2083 0.0417 0.0417 (37) 
0.0417 0.0417 —0.2083 0.0417 

0.0417 0.0417 0.0417 -0.2083 


Its associated quadratic Lyapunov function, V(x) based in (34) is given by: 


V(x) =- 0.2083x? +0.0834.x,x, + 0.0834-x,x, +0.0834.x,x, —0.2083.x; +0.0834.x,x, +0.0834x,x, 
+0.2083x; +0.0834.x,x, +0.2083.x; 
The eigenvalues of the matrix, P are given by: 
P(A) = [-0.2500, —0.2500, —02500, —0.0833] (38) 


2.5. Matrix exponential function analysis approach 
For any n x n square matrix, A that can be diagonalized as: 


A=SDS"1 (39) 


where matrix, S consist of eigenvectors of A, and D is a diagonal matrix with the eigenvalues of A along the 
diagonal, t is a real or complex variable; then its matrix exponential is given by: 


ease os (40) 


The (40) is expanded as follows: 


io) i 2 3 
gat = Čini aD aD L... (41) 


=1+4At+ + 
2! 3! 


Which is an imitation of the power series: 


= i (x)! _ x2 x3 
e~= ia ee tap (42) 
where, A, A?, A3, ...... are all an n x n matrix. A Laplace transform, L is an operator which takes a function 


F(t) as its input and produces f(s) as its input. The inverse Laplace Transform L~! takes f(s) as input and 
produces F(t) as output. Therefore, 
o (at)! - 
aeria siR - AI} (43) 
and 
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p(t) = LH{[sl, — ATt} = e“ (44) 
PE) = L*(sI — Ay) =I + tat 2 4 (45) 
(sf - A)" = (1/5) -4/51 = 1 44444 (46) 


where, (s7 — A) `! is known as the resolvent of A, L and L” is the Laplace transform and inverse Laplace 
transform respectively, and I, is n x n unit matrix. 

In (43) and (44) are applicable to time-varying linear systems and depends on the initial as well as 
the final time and not just the difference between them All linear systems that are asymptotically stable are 
also exponentially stable, [33], [34]. Applying in (39), (40), (43) and (44) respectively to the transition matrix 
given by (22): 


-1-1 -1 -1 
Z 1 0 0 1 (47) 
0 1 0 1 
0 0 1 1 
1 3 1 1 
4 4 4 4 
[i pa 3s. od (48) 
go | 4 4 4 4 
1 1 1 3 
4 4 4 4 
es tel te 
4 4 4 4 
2 0 0 0 
p= 0 2 0 0 (49) 
0 0 2 0 
00 0 6 
Substituting in (47)-(49) into (40) given by: 
-1 -1 -1 -1) fe” 0 0 0 1l 3 1 1 
pes 10 0 1 0 & 0 0 E s ; : 68) 
0 1 0 1| j0 0 &# 4 4 4 4 
00 1 1 00 0 & l 1 81 3 
4 4 4 4 
La Wey sed ee 
4 4 4 4 
3e” +e" e” pe” —e”pe" e” pe“ 
4 4 4 4 
—e” +e“ 3e” +e“ e“ e e“ — e” 5D 
gi 4 4 4 4 
e” — e” e“ Ze” e” +3e” e“ en 
4 4 4 4 
—e” +e“ onga e” — e” 3e” +e“ 
4 4 4 4 


A quick check to determine the accuracy of the obtained exponential matrix equation is to set t = 0 in the 
right hand side of (51) to obtain a unit or identity matrix, I given by (52) [35]: 
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100 0 
A) _] 9 10 0 (52) 
001 0 
000 1 


2.6. Hurwitz stability criterion approach 

Hurwitz criterion is an approach used to determine the stability of a system by using the coefficients 
of the characteristic polynomial equation without counting the roots involved. Given the characteristic 
polynomial [36]-[38]: 


p(A) = anA” + dyad t+. eee +a,A + ao (53) 


Where Qo,..-..+-0ee- ,An E R, Qao, An > 0. 
The characteristic polynomial, p(2) with real coefficients is stable or a Hurwitz polynomial if all its 


zeros have negative real parts. Hurwitz matrix, H for (53) is given by (54) [39], [40]: 


ip UL ie OS peer a 0 | 
a3 a9 ay Psst 0 (54) 
a a a One PA 0 
H, = 5 4 3 2 
ay a6 (E ee 0 
|f 2n-1 92n-2  %2n-3 | "Fn | 


The stability is determined using the sub-determinants of the characteristic polynomial, p(/) as follows [22], [38]: 


dD, = ay > 0 (55) 
aa 
D, = a = 1a, — agaz > 0 (56) 


The characteristic polynomial equation of the WT power generation system which is given by (29). 
Comparing in (29) with (53), the coefficients are: 


ao = 48, a, = —80, a, = 48, ag = —12,a, = 1 


The sub-determinants are therefore, given by: 


D, =a, = —80 (57) 
D, a = ay,az = aođ3 
es (58) 
D = | * °| = (—80)(48) — (48)(—12) = —3264 
a3Q2 
aa, 0 
D, = |@302@,| = azD, = (—12)(—3264) = 39,168 (59) 
00a, 


3. RESULTS AND DISCUSSIONS 
For a Lyapunov function in variables x1, x2,...... „Xn Which can be put in matrix form V(x) = x7 Px, 
where P is a real symmetric matrix; necessary and sufficient conditions for V(x) to be positive-definite are 
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provided by Sylvester’s criterion, which states that all the principal minors of P of order 1, 2, ....... „n must be 
positive. That is: 
PisPi2 P11P12P13 
P11 > 9, 3 > 0, |P12P22P23| > 0, -0000ean ,|P|>0 
TERSA P13P23P33 
Determining the principal minors based on Sylvester’s criterion, p;, = —0.2083. Thus, p4, < 0 does not 
meet the criteria. Pi1P12) _ 0.04165. Thus, P11P12 > Owhich meets the criteria. 
12P22 12P22 
P11P12P13 P11P12P13 
P12P22P23) = —0.00781. Thus, |P12P22P23| < 0 does not meet the criteria. 
P13P23P33 P13P23P33 
P11P12P13P14 P11P12P13P14 
P12P22P23P24) _ 0.0013. Thus, P12P22P23P24 Smeets: 
P13P23P33P23 P13P23P33P23 
P14P24P34P44 P14P24P34P44 


Thus, by Sylvester’s criterion, all the principal minors of P are not positive-definite and for that reason, the 
large-scale WT power system being considered is therefore, not asymptotically stable. 

For a system to be Hurwitz stable, the determinants must meet the following criteria: D7 > 0; D2 > 
0; and D3 > 0. The results obtained in the Hurwitz analysis are: D} = —80, D, = —3264 and D} = 39,168. 
From the results, only D3 meets Hurwitz criteria which is an indication that the large-scale WT power system 
under consideration is also not Hurwitz stable. In the case of the matrix exponential analysis, the large-scale 
wind power system is stable if e^t > 0 as t > 0. All trajectories of e^t must converge to zero (0) as t > 0. 


Thus, the obtained matrix A of the large-scale wind power system is stable if and only if all eigenvalues of 
matrix A have negative real part: RA; < 0,i = 1,2,n. As the lime“! + 0, is an indication that the large-scale 


WT power system as well is not asymptotically stable which means mode response decay cannot be 
guaranteed for large, t. The results obtained through this modelling approach and analysis is less time 
consuming as compared to computer modelling and simulations to observe the results. 


4. CONCLUSION 

Different approaches namely Lyapunov, exponential matrix and Hurwitz were used to determine the 
stability of the large-scale renewable energy power grid. All the three approaches confirmed that the 
stochastic model is unstable. The assumptions made during the modelling stage can be revised to improve 
upon the model’s stability; or a different modelling and analysis approach can be employed to ascertain the 
stability of such a power system. Ascertaining the actual behavior pattern will ensure that the appropriate 
control systems are designed for the stabilization of the WT power system based on the model. The primary 
contribution of this research is the development of a stochastic model for entirely renewable power grid 
which is envisaged to be implemented in the near future as result of dwindling fossil fuel reserves and 
environmental concerns. The stability of the grid was analyzed using three different approaches. All the three 
approaches used confirmed the hypothesis that the proposed large-scale renewable power grid might not be 
stable based on the modelling approach. It is noteworthy that this research did not factor battery storage 
system into the development of the model. 
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